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DEALING WITH UNCERTAINTIES IN INITIAL ORBIT 

DETERMINATION 

Roberto Armellin* Pierluigi Di Liziaj and Renato Zanetti* 


A method to deal with uncertainties in initial orbit determination (IOD) is pre- 
sented. This is based on the use of Taylor differential algebra (DA) to nonlinearly 
map the observation uncertainties from the observation space to the state space. 
When a minimum set of observations is available DA is used to expand the solu- 
tion of the IOD problem in Taylor series with respect to measurement errors. When 
more observations are available high order inversion tools are exploited to obtain 
full state pseudo-observations at a common epoch. The mean and covariance of 
these pseudo-observations are nonlinearly computed by evaluating the expectation 
of high order Taylor polynomials. Finally, a linear scheme is employed to update 
the current knowledge of the orbit. Angles-only observations are considered and 
simplified Keplerian dynamics adopted to ease the explanation. Three test cases of 
orbit determination of artificial satellites in different orbital regimes are presented 
to discuss the feature and performances of the proposed methodology. 


INTRODUCTION 

Orbit determination is typically divided into two phases. When the number of observation is 
equal to the number of unknowns a nonlinear system of equations need to be solved. This problem 
is known as initial (or preliminary) orbit determination (IOD). When more observations are available 
accurate orbit determination can be performed. IOD typically delivers a single solution (or a limited 
number of solutions) that exactly produces the available observations. In addition, in IOD simplified 
dynamical models are often used (e.g. Keplerian motion) and measurement errors are not taken 
into account (the problem is deterministic). When more observations are available the approach 
becomes stochastic, because the additional observations include noise. This problem is usually set 
as an optimization one, in which the (optimal) solution is the one that minimizes the observation 
residuals. The solution is obtained via batch estimation, e.g. weighted nonlinear least squares, or a 
sequential estimation, e.g. extended Kalman Filtering . 1 

In this paper we focus our attention on the orbit determination of resident space objects (RSO) 
observed on a single passage with optical sensors. Thus, the problem is the one of an angles- 
only orbit determination. In order to determine the orbit an IOD problem is solved followed by a 
procedure to update the initial solution based on the additional observations. 

Angles-only IOD an old problem. Gauss ’ 2 and Laplace’s 3 methods are commonly used to de- 
termine a Keplerian orbit that fits with three astrometric observations. These methods have been 
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revisited and analyzed by a large number of authors (e.g. 4,5, 6) and new ones introduced more 
recently. The Double r-iteration technique of Escobal 7 and the approach of Gooding 8 arc two ex- 
amples of angles-only methods introduced for the IOD of RSO. 

In 2012 Armellin at al. 9 proposed a IOD solver based on the solution of a Lambert’s problem 
(between the second and the third observations) and a Kepler’s problem between the first and second 
observation. The method iterates on the slant ranges at the second and third observations in order to 
drive to zero the observational defects at the first observation. The iterations were carried out with 
a high-order extension of Newton’s method enabled by differential algebra (DA). In addition, high 
order Taylor expansions were exploited to nonlinearly map the uncertainties from the observation 
space to the state space. 

In this work a modified version of the method is proposed, in which all the three slant ranges 
are the problem unknowns. The approach is based on the solution of two Lambert’s problems and 
using the continuity of the velocity vector at the central observation as constraint. The method 
has no restrictions on the geometry of the observations and it can deal with both short and long 
gaps. As in the previous work the solution is obtained with a high-order Newton’s iteration scheme 
enabled by DA. This approach allows the algorithm to both convergence in few iterations and map 
uncertainties form the observation space to the state space. Thus, already the initial orbit is provided 
with statistical information. 

When multiple observations on the same passage arc available the IOD solution is updated. In- 
stead of adopting a classical least squares approach (which employs the linearization of the dynam- 
ics and of the measurement functions 10 ) high order inversion tools available in DA arc exploited to 
nonlinearly map group of observations to the state space at a common epoch, thus producing full 
state pseudo-observations. The mean and covariance of these pseudo-observations arc nonlinearly 
computed by evaluating the expectation of the related high order Taylor polynomials. Linally, a lin- 
eal - updating scheme is utilized to update the current knowledge of the state mean and covariance. 

The paper is organized as follows. A brief introduction on the DA tools used for the implemen- 
tation of the algorithm is given first. This covers the methods to expansion the solution of ordinary 
differential equations (ODE), compute the expansion of the solution of implicit parametric equa- 
tions, and the algorithm to map statistics through nonlinear transformations. The following sections 
describes the main algorithms developed in this work, i.e. the angles-only IOD solver and the up- 
dating scheme. Simulated observational scenarios for a Geosynchronous Transfer Orbit (GTO), a 
Geosynchronous Orbit (GEO) and a Molniya are used to analyzed the performances of the imple- 
mented methods. Some final remarks conclude the paper. 

DIFFERENTIAL ALGEBRA TOOLS 

DA supplies the tools to compute the derivatives of functions within a computer environment. 1 1 
More specifically, by substituting the classical implementation of real algebra with the implemen- 
tation of a new algebra of Taylor polynomials, any function / of v variables is expanded into its 
Taylor polynomial up to an arbitrary order n with limited computational effort. In addition to basic 
algebraic operations, operations for differentiation and integration can be easily introduced in the 
algebra, thusly finalizing the definition of the differential algebra structure of DA. 12- 13 Similarly to 
algorithms for floating point arithmetic, also in DA various algorithms were introduced, including 
methods to perform composition of functions, to invert them, to solve nonlinear systems explicitly, 
and to treat common elementary functions. 14 The differential algebra used for the computations in 
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this work was implemented in the software COSY INFINITY. 15 The reader may refer to Di Lizia et 
al. 16 for the DA notation adopted throughout the paper. 

High-order expansion of the solution of ODE 

An important application of DA is the automatic high order expansion of the solution of an ODE 
in terms of the initial conditions. 14 - 16 This can be achieved by replacing the operations in a classical 
numerical integration scheme, including evaluation of the right hand side, by the corresponding 
DA operations. This way, starting from the DA representation of an initial condition xq, DA ODE 
integration allows the propagation of the Taylor expansion of the flow in xo forward in time, up to 
any final time tf. Any explicit ODE integration scheme can be rewritten as a DA integration scheme 
in a straight-forward way. For the numerical integrations presented in this paper, a DA version of a 
7/8 Dormand-Prince (8-th order solution for propagation, 7-th order solution for step size control) 
Runge-Kutta scheme is used. The main advantage of the DA-based approach is that there is no need 
to write and integrate variational equations in order to obtain high order expansions of the flow. It is 
therefore independent of the par ticular right hand side of the ODE and the method is quite efficient 
in terms of computational cost. 

Expansion of the solution of parametric implicit equations 

Well-established numerical techniques (e.g., Newton’s method) exist, which can effectively iden- 
tify the solution of a classical implicit equation 

f(x) = 0 (1) 

with / : — »• W . Suppose an explicit dependence on a vector of parameters p can be highlighted 

in the previous vector function /, which leads to the parametric implicit equation 

f(x,p) = 0. (2) 

Suppose the previous equation is to be solved, whose solution is represented by the function x(p) 
returning the value of x solving (2) for any value of p. Thus, the dependence of the solution of the 
implicit equation on p is of interest. DA techniques can effectively handle the previous problem 
by identifying the function x(p) in terms of its Taylor expansion with respect to p. This result is 
achieved by applying partial inversion techniques as detailed in 16. 

The final result is 

[x] = x + M x (6p), (3) 

which is the fe-th order Taylor expansion of the solution of the implicit equation. For every value 
of Sp, the approximate solution of fix. p) = 0 can be easily computed by evaluating the Taylor 
polynomial (3). Apparently, the solution obtained by means of Map (3) is a Taylor approximation 
of the exact solution of Eq. (2). The accuracy of the approximation depends on both the order of the 
Taylor expansion and the displacement dp from the reference value of the parameter. 

Nonlinear mapping of the estimate statistics 

Consider random variable x e ’ft" with probability density function p(x) and a second random 
variable y 6 3 ft m related to x through the nonlinear transformation 

y = /(*)• (4) 
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The problem is to calculate a consistent estimate of the main cumulants of the transformed probabil- 
ity density function p{y). Since / is a generic nonlinear function this formulation includes a wide 
range of problems involving uncertainty propagation (uncertainty propagation through nonlinear 
dynamics, uncertainty propagation through nonlinear coordinate transformations, etc.). 

The Taylor expansion of y with respect to deviations Sx can be obtained automatically by initial- 
izing the independent variable as a DA variable and evaluating (4) in DA framework. This procedure 
delivers 

[y] = /([as]) = y + My(Sx) = c Pl ,., Pn ■ Sx Pl ■ ■ ■ 5x Pn , (5) 

Pl-t \-Pn<k 

where in this expression y is the zeroth order term of the expansion map, and c pi ... Pn arc the Taylor 
coefficients of the resulting Taylor polynomial 

1 gpi+-+Pn f 

C pi--Pn pi\---p n ! dx p ' ■ ■ ■ dx P a ’ ^ 

The evaluation of (5) for a selected value of 5x supplies the k - th order Taylor approximation of y 
corresponding to the displaced independent variable. Of course, the accuracy of the expansion map 
is function of the expansion order and can be controlled by tuning it. 

The Taylor series in the form (5) can be used to efficiently compute the propagated statistics . 17 ' 18 
The method consists in analytically describing the statistics of the solution by computing the /-th 
moment of the transformed pdf using a proper form of the l - th power of the solution Map (5). 

For a generic scalar random variable x with pdf p(x) the first four moments can be written as 


/ 1 = E{x} 

P = E{(x- y) 2 } 

, E{(x-y) 3 } 

7 = ^3 

E{(*%) 4 } 


where ji is the mean value, P is the covariance, 7 and k are the skewness and the kurtosis, respec- 
tively , 19 and the expectation value of x is defined as 



xp(x)dx. 


( 8 ) 


The moments of the transformed pdf in (4) can be computed by applying the multivariate form 
of Eq. (7) to the Taylor expansion (5). The result for the first two moments becomes 

Py t = £{[?/,:]} = c i,Pi-Pn E { SxP i ■ ■ ■ 6x n n } 

Pl-t KPn<fc 

P y iVj = E {([yi] - Pi)([yj\ - Pj)} = C i,Pl-Pn C j,qi-qn E { Sx l 1+qi ■ ■ ■ Sx Pn+qn }, 

pi -I h Pn<k, 

qi-\ h qn<k 

(9) 

where c,;, pi ... Pn are the Taylor coefficients of the Taylor polynomial describing the i-th component 
of [y]. Note that in the covai'iance matrix formula the coefficients c PPl ... Pn and Cj m ... q . n are updated 



to include the subtraction of the mean. The coefficients of the higher order moments arc computed 
by implementing the required operations (e.g. Cy,\ — p i )(\y. J ] — fi :J ) for the second order moment) 
on Taylor polynomials in the DA framework. The expectation values on the right side of Eq. (9) 
arc function of p(x). It follows that if the initial distribution is known, all of the moments of the 
transformed pdf p(y) can be calculated. The number of monomials for which it is necessary to 
compute the expectation increases with the order of the Taylor expansion and, of course, with the 
order of the moment we want to compute. Note that, at this time, no hypothesis on the initial pdf has 
been made. Thus, the method can be applied independently of the considered variable distribution. 

We now consider the case in which x is a Gaussian random variable (GRV), x A f(li,P), 
in which p, is the mean vector and P the covariance matrix. An important property of Gaussian 
distributions is that the statistics of a GRV can be completely described by the first two moments. In 
case of zero mean, the expression for computing higher-order moments in terms of the covariance 
matrix is due to Isserlis. 20 In physics literature, Isserlis’s formula is known as the Wick’s formula. 


Let si to s n be nonnegative integers, and s = si + S2 + • • • + s n . Then the Wick’s formula 
suggests that 

/ 

0, if s is odd 


£T f /y 5 ! <y. S 2 r r S n\ 

^ l x l x 2 ' • • x n J ~ 


Haf(P), if s is even 


(10) 


where Haf(P) is the hafnian of P = (n r j), which is defined as 


Haf(P) = Yl a P2i-i,P2i , (U) 

perL i=1 

and f|.s is the set of all permutations p of {1,2 ,..., s} satisfying the property p\ < pz < Pb < 
... < p s -i and pi < p 2 , P3 < P4, ■ ■ ■ ,p s - 1 < Ps- 21 

We observe that the expectation value terms of Eq. (9) can be computed using Eq. (10), and the 
resulting moments can be used to describe the transformed pdf. 


DA-BASED ANGLES-ONLY IOD 

In the classical angles-only IOD problem three optical observations at epoch i,, with / = 1 .... . 3 
are available. The observations consist in three couples of right ascension and declination angles, 
(aii,di). These observations provide us with three inertial light of sights p,. i.e. the unit vectors 
pointing from the observer (on the Earth's surface) to the observed object. 

Assume to have first guess values of the slant ranges pi or equivalently for the orbit radii r t (e.g. 
from the solution Gauss’ 8-th degree polynomial). We present a high order iterative procedure with 
the following objectives: a) refine the values of pi assuming Keplerian dynamics, and b) express the 
functional dependence of the solution of the IOD problem with respect to observation uncertainties 
in terms of a high-order Taylor polynomials. 

We start by initializing the observations as DA variables: 

[a] = a + 5a 
[5]=<5 + M, 

in which we have grouped the observations in two homogeneous vectors, a = (01,02,0:3) and 
<5 = (5i, 62 , A, ) , and 5a and 65 accounts for measurement uncertainties. The line of sight vectors 
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at 1 1 , 1 2 and t : > become 


[Pi\ = Pi+Mp^SauSSi) 

[p2[ = P2 + -Mp 2 (5a 2 , M 2 ) (13) 

[As] = A3 + -Mp 3 ( 6 a 3 , 663 ), 

where M.p. is an arbitrary order Taylor polynomial that describes the effect of an observation un- 
certainty on the line of sight. 

Similarly, we initialize DA variables on the topocentric distances at t\, t - 2 and £3 


M 1 = Ai + $Pi 

[A 2] 1 =A 2 +^A 2 (14) 

[A3] 1 ” = A 3 " + 5 A3, 


or in more compact form 

[ P r =p i ~+5 P , as) 

where the superscript 1 ” indicates the first step of the iterative procedure, and p\ , p\ , and p\ 
are the guess values for the slant ranges. 

The spacecraft position vectors can be written (by summing the known observer’s locations) as 


[r 1 ] = tt + M ri (5a 1 ,55 1 ,5pi) 

[r 2 ] = r 2 + M r 2 (5a 2 ,55 2 ,5p 2 ) (16) 

[r 3 ] = r 3 + M r3 (5a 3 ,65 3 ,Sp 3 ). 


A DA-based Lambert’s problem 22 can be solved between with [r*i] and [r 2 ], and between [r 2 ] 
and [r 3 ]. Using the DA-implementation of Lambert’s problem we obtain two Taylor polynomial 
approximations for the velocity vector at t 2 

[ v ~] = v - + M v -(5ai,55i,5a 2 ,55 2 ,5pi,5p 2 ) 

Hi = v% + M v +( 6 a 2 ,S 6 2 , 5 a 3 , 563 , 6 p 2 , 5 p 3 ) 1 

Note that the above expressions of the velocity vector arc different for two reasons. First, the 
starting values of the slant ranges arc not the solution of the IOD problem; secondly, they have 
different functional dependence on the observation angles. The goal is thus a) to find the values of 
the slant ranges such that the velocity vector is continuos at the midpoint, i.e., we want to find the 
exact values of pi, p 2 , and p 3 , and b) to approximate the spacecraft state at t 2 as a Taylor polynomial 
in the observation uncertainties. We start by defining the Taylor map of the defects 

[AF 2 ] = Hi ~ H] = ^2 + M&v 2 (5u, 55 , 5p). (18) 

Note that for the exact values of p \ , p 2 and p 3 the constant paid of maps (18) would be zero. We 
now need to find the variations 5p necessary to cancel out these constants and to express r 2 and v 2 
as Taylor polynomials in 5a and 65 only. The first step is to work with an origin preserving map 

[At> 2 ] = [Aha] - Av 2 = MAv 2 (<> a , 5p) (19) 
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and to build an augmented Taylor polynomial by adding identities in observation deltas 


1 

> 

£ 

to 

1 


-Ad a V 2 


6 a 

6 a 

= 



66 

5<5 


Is 


. V . 


(20) 


This polynomial map can be inverted using ad-hoc algorithms implemented in COS Y-Infinity, yield- 
ing 


6 a 


-Ad Av2 

- 1 

Av-2 

66 

= 

I a 


6 a 

_ 6 P . 


Is 


66 


Extracting the three last lines we obtain 


( 21 ) 


[ 5P 


; m p ; 


At>2 

5 a 

66 


( 22 ) 


We now evaluate the map ( 22 ) in = —Av2, obtaining 


[p\ 1+ = p 1+ + M p ( 5 a, 66) 


( 23 ) 


where the subscript 1 + indicates the Taylor polynomial of the corrections of the topocentric dis- 
tances to be applied at the end of the first iteration. This last step is the high-order counterpart of 
classical Newton’s method. 


The second iteration starts with the Taylor polynomials of the topocentric distances given by 

[/ P ] 2 =[p\ 1 +[p] 1+ +6p = p 2 + M p ( 5 a, 56 , 6 p) ( 24 ) 

where now the explicit dependence on the entire set of observables appears. Thus, from the second 
iteration, the Taylor polynomials ( 16 )— ( 17 ) depend on all (6a, 66, 6 p). The iterative procedure 
ends when the values of A62 arc smaller than a prescribed tolerance. The Taylor polynomials of 
the topocentric distances at the last iteration k arc 

[p\ = [p\ k ~ + [p\ k+ = P + M p (6a, 66) ( 25 ) 


Using these expressions the spacecraft position and velocity vectors at t-2 assume the form 



H 

= r 2 + M r2 (6a,66) 

( 26 ) 


M 

= v 2 + M V2 (6a, 66). 

or more compactly 

[*2] 

= x 2 + M X2 (6a, 66), 

( 27 ) 

where X2 = (r 2,1)2)- 





As a result of the iterative procedure, r-2 and v-2 exactly satisfy (in the two-body model) the nom- 
inal observation set (a, 6). Furthermore, for any displaced value of the observables, the solution of 
the preliminary determination problem is computed by evaluating the polynomial ( 26 ) in the corre- 
sponding values of (da, 66). Map ( 27 ) is an arbitrary order Taylor polynomial in 6 a and 66 , which 
maps the uncertainties from the observable space to the spacecraft state phase. In particular, using 
the approach described in Section “Nonlinear mapping of the estimate statistics” we can compute 
the statistical moments of x, given the statistics of the measurements. 
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DA-INVERSION IOD 


When more than three optical observations arc available the solution (reference state and associ- 
ated statistics) of the IOD problem needs to be updated to include the additional information. This 
is carried out through a high-order filtering technique based on nonlinear mapping of statistics and 
lineal - update scheme, in which only the pdf of the measurements is constrained to be Gaussian. 

The optimal linear estimate of a state x based on a measurement y is given by 

* = /*» + p xyPyl {y~ p v ) ( 28 ) 

where /i,,. is the state mean, P xy is the joint covariance of the state and the measurement, and 
P yy is the covariance of the measurement. For a general non-linear measurement with additive 
noise y = h(x) + r/, calculating fi y and the covariance matrices requires full knowledge of the 
distribution of the state. This requirement has two consequences: first it means that the state and 
its uncertainty need to be propagated forward to the measurement time, and second that statistics 
of the measurement need to be calculated through a nonlinear transformation of the current state. 
In this work we propose addressing this issue in a different way. The state is always estimated at a 
fixed epoch time, and the nonlinear map to transport it to any other epoch is calculated with the DA 
framework. Instead of working with y as a function of x, a full pseudo-measurement of the state 
is generated from y; the inverse of the non-linear map from the state to the measurement is readily 
available from COSY-Infinity. The advantage of this approach is that only the distribution of the 
measurement noise is assumed Gaussian while the distribution of the state is left unconstrained. 

Consider a time span [to, tf] and let x k be the state variable at some time t/ c G [to, t/]. Consider 
also a set of N measurements y i given at times t, <£ [to, f/] with i = 1, . . . , N. Given the current 
estimate of the state xT and the related error statistics, we can always define the estimated state 
as a DA variable and compute the predicted measurement at t, in the DA framework. The relation 
between state and measurement is a nonlinear map that accounts for the forward propagation of the 
initial condition and the measurement function. Under proper conditions this relation can be inverted 
to map the observation space at into the state space at tk- The main cumulants of the resulting 
map can be computed as described in the previous section, with the assumption that the statistics 
of the measurement errors is Gaussian. The computed mean and covariance are exploited to update 
the knowledge of x k using a linear update scheme. This can be done for groups of measurements 
for which the dimension of measurement vector y, is equal to the dimension of the state vector, and 
the map is invertible. 

The resulting method can be made recursive and summarized as follows. From the IOD algorithm 
we start from an initial value of the state estimate and covariance, x y = y,~ k and P XkXk (* n general 
tk = t- 2 , the epoch of the central observation in the IOD problem.) Define the current estimate at 
time of interest tk as a DA variable; i.e.. 


[x k \ = x k + 5x k . (29) 

and propagate it to time t, when a measurement becomes available. The result assumes the form of 
the following high-order Taylor expansion map 

[x.i] = Xi + M Xi (6x k ). (30) 

Note that the constant part of this map, i.e. x t , it is not the predicted mean at time t, due to the 
nonlinearities of the dynamics (the relation x t = y,, r holds true only if the state transition matrix is 




Figure 1: Sketch of the Taylor maps involved in the construction of the DA-base map inversion 
nonlinear filter. 


used). Then, use the measurement equation to compute 

[Vil = &([*»]) = Vi + M y .(5xk), (31) 

where h represents the measurement function. Figure 1(a) can be used by the reader to better 
understand the meaning of Maps (30)— (3 1). 

The next step consists in defining an origin preserving map 

= [Vi] ~Vi = Msy.{8x k ). (32) 

This polynomial map can be inverted if two conditions are satisfied: the map must be square and all 
the measurements must be independent. If these requirements arc satisfied, we can invert Map (32) 
using algorithms implemented in COSY-Infinity, obtaining 

5x k = Msx k (8yi)- (33) 

We now substitute in Map (29) the expression of 5x k from (33), yielding 

[x k \ = x~ + M Xk (8yi). (34) 

This map now represents the pseudo-measurement of state x k based on the observation y v so it is 
renamed as 

l z k) = + -A'Wfo/j). (35) 

By construction the constant paid of Eq. (35) is equal to the state estimate at step k, i.e. x^, but 
its statistical moments arc different to those of x k , due to the nonlinear contribution of M. Zk (6y i ) 
(as highlighted in Fig. 1(b)). We can now apply Eq. (9) to Taylor expansion (35) to compute the 
statistics of the random variable z k and, in particular, the first two moments y, Zj and P ZkZk . The 
computed mean can be treated as the “predicted measure” of the state at time t k , with measurement 
error defined by P ZkZk . Thus, we can update the initial estimate and error covariance, using the 
least squares method. This can be done using the Kalman filter update equations that, applied to the 
current problem, read 


K=p- kXk (P~ kXk +P ZkZk )-\ (36) 

H +K(z k - y Zk ) , (37) 

Pt kXk ={I~K) P~ kXk (I - Kf + KP ZkZk K T , (38) 
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where x £ is the updated mean at time tf~ and Px k x k the related updated covariance matrix. When 
another measurement becomes available, we can define the state at time ty- as a new DA variable, 
centered in the new estimate i it, and iterate the process. Note that zy. is the true state-measurement 
at t, mapped to time A, which is readily available by evaluating Map (35) for 5y i = ij, — y,. 

We said that Map (32) must be square in order to be invertible. It follows that if the measurement 
vector has smaller dimension than the state vector, after the first measurement is received we can not 
proceed with the update, but we have to wait for additional measurements (i.e. in the optical case 
three observations arc needed). When the number of scalar measurements equals the dimension of 
the state variable, we can define an augmented measurement vector that can be used to build Maps 
(31) and (32). 

Once the final estimate of the state at time ty- is obtained, the statistics of the solution can be 
computed at any time via propagation and DA-based expectation evaluation. 

TEST CASES 

The algorithms for IOD arc run considering single -pass optical observations of three objects as 
listed in Table 1 . 


Table 1: Test cases: orbital parameters 


Test Case 

A 

B 

C 

Orbit type 

GEO 

GTO 

Molniya 

NORAD ID 

26824 

23238 

40296 

Epoch 

JED 

2457163.2824 

2457167.1008 

2457165.0708 

a 

km 

42143.781 

24628.972 

26569.833 

e 

- 

0.000226 

0.699849 

0.723221 

i 

deg 

0.0356 

3.962 

62.794 

n 

deg 

26.278 

315.676 

344.538 

UJ 

deg 

42.052 

240.885 

271.348 

M 

deg 

72.455 

13.735 

347.726 


The observations arc all simulated from Teide Observatory, Tenerife, Canary Islands, Spain (ob- 
servation code 954). The simulation windows are summarized in Table 2. For all the cases 15 
equally spaced optical observations arc simulated within the observation window. The spacecraft is 
considered observable when its elevation is above 10 deg, is in sunlight, and the Sun has an eleva- 
tion lower than -7 deg. As a result, different observation gaps are considered, ranging from 522 s 
for the GTO case to 2160 s for the GEO case. The GTO object is observed before the apogee for an 
arc length of approximately 20.7 deg. The average separation between observations is 1.5 deg, with 
maximum and minimum values of 1.9 and 1.3 deg, respectively. The Molniya object is observed 
before the apogee on an arc length of 13.4 deg. In this case the mean, maximum, and minimum 
observation separations are 1, 1.1, and 0.8 deg. Finally, for the GEO case the observed arc has a 
length of 127.4 deg (with uniformly spaced observations). 

For all the cases the central observations, i.e. observation ID 7, 8, and 9, are used for the IOD; 
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thus, ccg = (fg, Eg) and P~ 8 Xg are the output of the IOD problem. The remaining observations arc 
used for the update of rcg and P~ 8 . Finally, pertaining to the accuracies, we consider Gaussian 

measurement noises with standard deviation of 0.5 arcsec. 


Table 2: Test cases: observation windows 


Test Case Observation Window At cr a g 



yr 

mo 

day 0 

day/ 

hr 0 

hr / 

hr 

arcsec 

A 

2015 

MAY 

22 

23 

21.000 

05.400 

0.600 

0.5 

B 

2015 

JUN 

02 

02 

03.550 

05.580 

0.145 

0.5 

C 

2015 

MAY 

22 

22 

20.600 

23.400 

0.200 

0.5 


All simulations are run on a MacBook Air with a 1.8 GHz Intel i5 CPU and 4 GB RAM. 


DA-based angles-only IOD 

The IOD algorithm is run 100 times for each of three test cases described in Tables 1 and 2. 
The observation geometries arc described in Figures 2(a), 2(c), and 2(e). For all the cases 6-th 
order computations are carried out. The DA-based IOD algorithm converges in all cases in, on 
average, three iterations (convergence is achieved when the euclidean norm of the velocity vector 
discontinuity at the central observation is less than 1 x 10~ 12 km/s). In all cases, the real solutions 
of the Gauss’ 8-th degree polynomial arc taken as first guesses for the unknown slant ranges. 

The result of the DA-based IOD algorithm is the Taylor polynomial [atg] (see Eq. (27)) which 
maps the observation uncertainties into uncertainties in the state space. This map is employed 
to compute the starting state estimate Xg and covariance P“ ga , g , evaluating the expectation of the 
monomials by assuming Gaussian statistics for measurement noise. Figures 2(f), 2(d), and 2(f) show 
the absolute value of the observation residuals associated to Xg (normalized by the observations 
standard deviation) at the different observation epochs and for all the 100 simulations. As expected 
the residuals arc minimal at the epochs of the IOD (i.e. ID 7, 8, and 9), whereas they steeply 
increase far from the central observations. In addition, note that Xg does not exactly satisfy the 
IOD, as it is acually the constant part of the associated Taylor polynomial, [reg], that does it (with 
an accuracy that depends on the threshold selected for algorithm convergence). The maximum 
differences between the constant paid of the map and the computed mean arc given in the first two 
columns of Table 3, where the contributions arc split in position and velocity components. It is 
apparent that the nonlinearities play a minor role for the test case A, and this is confirmed by the 
fact that the residuals arc minimal at observations 7, 8, and 9 for this test case (see Figure 2(f)). 

In all the cases the estimated covariance P~ ga , g is stretched along the line of sight directions as 
shown in the zoomed portions of Figures 2(a), 2(c), and 2(e). Higher nonlinearities affect test cases 
B and C, for which the uncertainty set is much more stretched. To quantify this, the maximum of 
the square root of the position and velocity eigenvalues (indicated with max <j rg and max a V8 ) arc 
reported in Table 3. 
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Table 3: IOD: uncertainty set description. 


Test Case 

max rg — f 8 1 1 

km 

max |«g — Fg | 
m/s 

max a rg 
km 

max a vs 
m/s 

A 

0.045 

0.003 

26.528 

1.976 

B 

7.579 

0.349 

340.993 

14.611 

C 

22.435 

1.312 

573.765 

30.675 


DA-based inversion IOD 

The results obtained by applying the updating scheme presented in Sec. “DA-inversion IOD” arc 
presented in this section. 100 simulations arc run for each of test cases and all the computations arc 
carried out at order 6, as for the DA-based IOD. 

As we arc considering 15 equally spaced optical observations, the maximum number of iterations 
(including the IOD using observations 7, 8, and 9) is 5. The updating scheme is stopped whenever 
the maximum number of iteration is reached or when the variation in the estimated state gets bigger 
than 5 times the maximum eigenvalues of the starting state covariance (this is considered as an 
anomaly in the updating scheme). 

For all the cases a set of 4 plots is presented. In the first one the difference between the current 
state estimate and the true state (indicated as 1 1 r 8 — 1 1 for position and 1 1 vg — \ | for velocity) is 

plotted as function of the iteration number. Mean, maximum and minimum values for the considered 
100 simulations are shown with different markers. In the second figure the maximum (over the 100 
simulations) of the maximum position and velocity eigenvalues of the estimated covariance matrix 
arc plotted as a function of the iteration number. Thus, the first two figures can be used to extract 
informations on state accuracy estimation and size of the estimated final uncertainty set. The third 
and fourth figures are about the observations residuals. More specifically, in the third figure the 
evolution of the mean residuals with the iteration number is highlighted using markers in gray scale 
(black markers for the last iteration); whereas in the fourth figure we plot the mean, maximum, and 
minimum values of the residuals (absolute value) at the the fifth iteration only. 

Figures 3, 4, and 5 show all a similar behaviour of the relevant quantities. The accuracy of the 
estimation improves with iteration number, and the size of the estimated state covariance reduces 
accordingly. The observation residuals decrease and become more homogeneous with the iteration 
number. More accurate predictions arc obtained for the Test Case A, thanks to both a longer ob- 
served arc and lower eccentricity of the orbit. In this case all the 100 simulations reach the 5-th 
iteration, with a mean final average estimation error of 0.164 km on position and 0.022 m/s on ve- 
locity. These errors increase to 3.353 km and 0.439 m/s for Test Case B, and to 8.520 km and 1.481 
m/s for the Test Case C. Note that the 96% of the simulations reach the fifth iteration for the for the 
GTO case, and this number further reduces to 90% for the Molniya orbit. 

Finally, in Figure 6 the results of 100 simulations using first order Taylor expansions are shown to 
highlight the effect of nonlinearities. It can be noticed that for the GEO case (Figure 6(a) and 6(b)) 
the updating algorithm is still convergent (although the average estimation error doubles with respect 
to 6-order expansion) as both the estimation errors and the residuals decrease with the iteration 
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number. This is not the case for both Test Case B and C, where the estimation errors and residuals 
decrease only up to the third iteration (i.e. when nine optical observations arc used). Thus, in 
these cases a linear approximation is not sufficiently accurate in mapping, to the central epoch, the 
observations taken at the boundary of the visibility windows. 

CONCLUSIONS 

In this paper the problem of dealing with observation uncertainties in IOD is addressed. A fully 
nonlinear method for IOD is implemented based on the high Taylor expansions delivered by DA 
computation. The method, based on the solution of two Lambert’s problems, delivers the solution 
of IOD problem and nonlinearly maps uncertainties from the observations space to the state space 
already when the minimum (three) number of optical observations arc considered. The algorithm 
converges for all the cases considered within, on average, three iterations. The average computa- 
tional time is 3.6 s when 6-th order computations arc carried out. 

A lineal - scheme for updating the state’s first two statistical moments is proposed when more 
optical observations are available in a single passage. This scheme is based on the generation 
of full state pseudo-observations at a common epoch, taking advantage of polynomial inversion 
tools available in DA. The required expectation are computed on high order Taylor polynomials, 
limiting the Gaussian assumption to the observation noises only. The updating schemes is shown 
to improve the accuracy of state estimation when short-dense observation arcs are available. The 
average computational time for the updating scheme is 1.91 s at order 6. 

In the present work simplified Keplerian dynamics are used. The algorithms can be easily ex- 
tended to arbitrary dynamics by using the DA-based tools for the Taylor expansion of the solution 
of ODEs (see 17 for details) and by replacing the Lambert’s solver with a DA-based algorithm for 
expanding the solution of two-point boundary values problems (as illustrated in 16). The authors 
plan to apply the algorithms to real observations including the case of short-dense radar observa- 
tions. 
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